Hybrid phase-space simulation method for interacting Bose fields 
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We introduce an approximate phase-space technique to simulate the quantum dynamics of inter- 
acting bosons. With the future goal of treating Bose-Einstein condensate systems, the method is 
designed for systems with a natural separation into highly occupied (condensed) modes and lightly 
occupied modes. The method self-consistently uses the Wigner representation to treat highly occu- 
pied modes and the positive-P representation for lightly occupied modes. In this method, truncation 
of higher-derivative terms from the Fokker-Planck equation is usually necessary. However, at least 
in the cases investigated here, the resulting systematic error, over a finite time, vanishes in the limit 
of large Wigner occupation numbers. We tested the method on a system of two interacting an- 
harmonic oscillators, with high and low occupations, respectively. The Hybrid method successfully 
predicted atomic quadratures to a useful simulation time 60 times longer than that of the positive-P 
method. The truncated Wigner method also performed well in this test. For the prediction of the 
correlation in a quantum nondemolition measurement scheme, for this same system, the Hybrid 
method gave excellent agreement with the exact result, while the truncated Wigner method showed 
a large systematic error. 

PACS numbers: 03.75.-b, 05.10.Gg, 02.50.Fz, 34.50.-s 
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I. INTRODUCTION 

The aim of this paper is to introduce a new, approx- 
imate, stochastic phase-space method and to test it on 
some simple problems with interacting Bose fields. A fu- 
ture goal of our research is to use the method to simulate 
the dynamics of interacting Bose-Einstein condensates 
(BECs). The method is, in fact, designed for BEC prob- 
lems, since it relies on the ability to make a meaningful 
separation of a multimode system into highly occupied 
(condensed) modes and lightly occupied modes. Hence 
our two-mode test cases will be constructed to have one 
highly occupied mode (N ^> 1) and one lightly occupied 
mode (N < 1). 

Previous authors have developed formalisms in which 
condensed atoms are treated in a different way to non- 
condensed atoms. Castin and Dum [lj studied the dy- 
namics of Bose-Einstein condensates at very low temper- 
atures using a Bogoliubov [2| approach, in which the bo- 
son field operator is written as a sum of condensate-mode 
terms and non-condensate-mode terms. Their treat- 
ment deals with number eigenstates rather than coher- 
ent states. They obtain results as an asymptotic expan- 
sion in the square root of the fraction of non-condensed 
atoms. Gardiner and Zoller in the third of their 
series on quantum kinetic theory, consider a stationary 
non-condensate band at a fixed temperature acting as a 
reservoir to the dynamic condensate modes. They de- 
rive a master equation for the condensate modes, using 
a number-conserving formalism. Dalton 0| calculates 
quantum correlation functions for boson field operators 
to use in the interpretation of double-well BEC inter- 
ferometry experiments. The approach is a phase-space 
method for a distribution functional, in which the Wigner 
representation is used for the condensed modes and the 
positive-P representation for the non-condensed modes. 



Our method will be seen to be substantially different from 
these three approaches. 

Besides in BEC evolution and collision problems, other 
typical cases where disparate occupation numbers exist 
would be in the quantum Brownian motion of a small 
number of massive particles inside a BEC, or in the col- 
lision of weak and strong coherent light pulses in a non- 
linear optical fibre. Hence we also consider these systems 
to be candidates for the Hybrid method. 

The foundations of this work are the stochastic phase- 
space methods developed to simulate the quantum dy- 
namics of systems with many degrees of freedom. In par- 
ticular we consider the Wigner-Moyal @, [f| approach, 
and the positive-P method [3, H]. We will see that both 
methods have wide applicability, but are ultimately lim- 
ited in the parameter regimes on which they can be used. 
The Wigner-Moyal method generally requires a trunca- 
tion to be able to map to a stochastic process. The result- 
ing approximate theory typically fails to give correct re- 
sults when significant numbers of modes with small mode 
occupation numbers are present 

The positive-P method is exact, but when applied to 
large multimode problems can often be used only for lim- 
ited simulation times before very large sampling error 
renders it unusable. The longest useful simulation times, 
for a given interaction strength, are for lightly occupied 
modes 

The new phase-space method to be introduced here 
is a combination of the Wigner and positive-P methods. 
In this Hybrid method, as we will call it, highly occu- 
pied modes are treated with the Wigner representation 
while lightly occupied modes use the positive-P represen- 
tation. A truncation of higher-order derivative terms is 
usually needed, but the resulting approximate method is 
expected to be valid (over finite times) to within correc- 
tions of the order of the reciprocal of the large occupation 
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numbers. 

The Wigner method is used in the regime where it is 
known to perform best and produces most simplification 
of the stochastic differential equations. The positive-P 
method is used on the modes that introduce most error 
in the truncated Wigner method. This latter choice is 
also designed to lengthen the useful simulation time. 

In this paper we will summarize the properties of the 
two representations, and discuss their successes and prob- 
lems, before actual construction of the Hybrid method. 
As a test case, we will apply the method to an exactly 
solvable problem: a system of two coupled anharmonic 
oscillators, one highly occupied, the other lightly occu- 
pied. The interaction preserves individual particle num- 
bers. 

At first we simply calculate the expectation values 
of quadratures and compare with the truncated Wigner 
method, the positive-P method and the exact solution. 
Then we investigate a higher-order correlation in the 
same system, one that would be observed in a quantum 
nondemolition measurement (QND) scheme. 



equation (|2.ip has a unique inverse, defining the Wigner 
function in terms of the density matrix: 

W(a) = J ^| e (-« Q *+«* Q ) Ti-(pe ( « at -«* a) ). (2.3) 



By manipulating equation (|2,2p . these basis operators 
can be written in the normally ordered Gaussian form of 
Corney and Drummond [12], 



A w {a) = 2 : e -*& -c?K&-<*) 



(2.4) 



where : f(a,a>) : indicates normal ordering. 

Now we may define the doubled Wigner representation 
with an expansion of the density matrix of the form 



p = J d a J d' ! a + W(a,a + )A w (a 1 a + ). (2.5) 

Here W(a, a + ) is a Wigner function defined on a doubled 
phase space and 



-2(a t -a+)(a-a) . 



(2.6) 



II. THE SINGLE AND DOUBLED WIGNER 
REPRESENTATIONS 

We consider a quantum many-body system of bosons. 
The relevant creation and annihilation operators are de- 



noted at, 



In the Wigner-Moyal approach, one com- 



plex phase space variable, a m , is used for each mode, m, 
of a system, and we call this a single phase space. In con- 
trast, a doubled phase space uses two complex variables, 
a m and a+, for each mode. We will find that using 
the Wigner and positive-P representations for different 
modes of the same system will generally require using 
a doubled phase space, although this can be avoided in 
certain cases. 

We begin by showing the definition and properties of 
the doubled Wigner representation. This is an extension 
of the familiar single phase space Wigner representation 
to a doubled phase space, and has been studied and ap- 
plied by Plimak et al Throughout this paper we will 
set h = 1. 

The single phase-space Wigner representation of the 
density matrix is given by 



P 



d 2 aW(a)A w (a). 



(2-1) 



This is an expansion of the density matrix on a basis 
of operators, the standard form we will use to compare 
all representations. Here W(a) is the Wigner function 
on phase space. Following Moyal [(| and Glauber and 
Cahill 0, 



A w (a) 
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are the operator basis elements, also defined on the dou- 
bled phase space. The new variable a + appears where a* 
had been, but in a stochastic simulation may take values 
different from the complex conjugate of a. 

From equation l|2.6p we can derive the operator corre- 
spondences for the doubled Wigner representation. The 
action of a creation or annihilation operator, multiplying 
the density matrix to the left or right, is equivalent to a 
linear differential operator acting on the Wigner function: 



ap 



pa 



1 d 



2da+ 



1 d 



W(a,a+) 



"'/' — | <> " 2 9a ' W(a,a + ) 



pa' 



2 oa 



(2.7) 



(2i 



(2.9) 



(2.10) 



We add a cautionary note. The derivation of equa- 
tions (|2.7H2.10p depends on the vanishing of boundary 
terms in an integration by parts. This problem is dis- 
cussed in Section IV. 

We note that a pure coherent state (with p = I7X7I) 
can be represented, in the doubled Wigner representa- 
tion, with the stochastic prescription 



1/ 

a = 7 + -(ni + in 2 ) 



(2.11) 



is an operator function on phase space, with trace unity. 
We also refer to this as the operator basis. We note that 



a + = 7* + -(ni - in 2 ) 



(2.12) 
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where m and n-i are independent real Gaussian random 
noises with unit standard deviations. 

In this symmetrically ordered representation, the for- 
mula for estimating symmetrically averaged products of 
creation and annihilation operators as stochastic averages 
over trajectories is: 

(a+ m a") sym = ((a +m a n )). (2.13) 

We will use the notation (()) throughout to signify a 
stochastic average over an ensemble of trajectories. 

We note that we will be exploiting the nonuniqueness 
of this representation and that of the positive-P repre- 
sentation: an infinity of different functions W(a, a + ) can 
give the same density matrix according to equation H2.5J1 • 
This feature of representations on doubled phase spaces 
will allow us, in the case of the Hybrid representation, 
to construct quasiprobabilities that are everywhere real 
and non-negative, obeying Fokker-Planck equations that 
allow mapping to a stochastic simulation. 



III. PROBLEMS WITH THE TRUNCATED 
WIGNER METHOD 

The (single phase space) Wigner representation has 
been widely used to study diverse physical problems, with 
great success [ijj, H, LLll ■ But a truncation of terms is 
necessary in most applications to allow a stochastic simu- 
lation. The truncated Wigner method is then not exact. 
Although the method gives good results in many cases, 
because the truncation of terms can be well justified if 
mode occupation numbers are large and simulation times 
are limited, the systematic errors can be significant if 
those conditions are not met. In addition, the estimation 
of a higher order moment (the expectation of a product 
of more than one field operator) will generally contain a 
larger systematic error than the estimation of the expec- 
tation value of a single field operator [161 ] . 

Even when no truncation is necessary, there is the 
problem of large sampling error in a truncated Wigner 
simulation. While an initial coherent state can be rep- 
resented by a positive-P distribution of zero width (see 
equation i|4.10p ). the same state will have a Wigner distri- 
bution with a finite width (equations l|2.1ip . (|2.12p ). For 
short times, the growing positive-P noise will not over- 
take the relatively constant Wigner noise. The result is 
greater sampling error in the Wigner simulation, requir- 
ing the calculation of far more trajectories to achieve the 
same precision. 

The investigation of Deuar and Drummond [§] into 
BEC scattering showed how the truncation problem pro- 
duces serious systematic errors in the simulation of a 
large number of interacting modes with many lightly oc- 
cupied modes. We will discuss these problems later in 
this section. 

Here we outline the reasons for truncation and the re- 
gion of validity of the approximation. 



From equation l|2.3p it may be seen that the Wigner 
function is always real, but it may take negative values 
for some density matrices. This would prevent us from 
mapping our quantum mechanics problem to a stochas- 
tic simulation, since the latter would require a positive 
semidefinite quasiprobability distribution. 

However, when we find the equation of motion for the 
Wigner function, the opportunity for an approximation 
procedure becomes apparent. This equation follows from 
the operator correspondences of the single Wigner rep- 
resentation (which can be obtained from equations <|2.^ 
I2.10p with the replacement a + — > a*) and the evolution 
equation for the density matrix 

^ = -i[H,p]. (3.1) 

We are going to restrict our attention to Hamiltonians, 
including multimode Hamiltonians, that include prod- 
ucts of creation and annihilation operators only up to 
quartic terms. This restriction will include the model of 
BECs with two-body s-wave scattering [T3|. The equa- 
tion for the evolution of a Wigner function under such a 
Hamiltonian will always take the general form 

dW d d 

~dT = -fo( A (<x)W(a)) - —(A*(a)W(a)) + T 3 . 

(3.2) 

Here T 3 is a term with three derivative operators, each 
either ^ or The key point to note is that for un- 
damped (unitary) time-evolution, there are never any 
second order (diffusion) terms, which is a consequence of 
the fact that the Wigner representation is symmetrically 
ordered. Also, fourth-order terms always cancel. These 
general results for quartic Hamiltonians, for the Wigner 
representation and for the positive-P representation, are 
summarized in Table 1. 





Drift Terms 


Diffusion Terms 


Third- Order Terms 


Wigner 


Yes 


No 


Yes 


Positive-P 


Yes 


Yes 


No 



Table 1: Terms in the Fokker-Planck equation for a quartic 
Hamiltonian, using the Wigner and positive-P 
representations. 



It is found that the third order terms, including for 
more general multimode problems, may be truncated and 
produce a systematic error in that is relatively small 
compared to the other terms, in the limit that the occu- 
pation numbers of the modes remain very much greater 
than unity. 

The motivation for this truncation is clear: equa- 
tion ll3.2tl then reduces to Liouville form, a special case 
of the Fokker-Planck equation in which only drift terms 
influence the evolution of the quasiprobability. If the 
initial density matrix for the problem is such that the 
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Wigner function is everywhere non-negative (and this is 
a common situation) then the function will remain non- 
negative for all times. A further mapping to a stochas- 
tic simulation becomes possible. The only noise in the 
simulation will come from the initial condition, since no 
second-order terms are present to cause diffusion. 

A small error in will produce a large error in W 
after a sufficiently long time, so this approximation pro- 
cedure can only be valid for a finite time. Over the rel- 
evant time-scales, the truncation is justified by a scaling 
argument. If the stochastic variable a is seen from the 
truncated equations of motion to remain of very large 
magnitude (\a\ ~ \/No ^> 1), then we define a scaled 
variable z — a/s/No and find that the third-order terms 
take the form 



t 3 ~ ±ddd((;w), 

ivo 



(3.3) 



where d is either S- or jtt and C is either z or z* . 

az az i ^ 

Deuar and Drummond [9] applied the truncated 
Wigner method to the large multimode problem of scat- 
tering BECs and found an ultraviolet divergence prob- 
lem: systematic errors that grow with the momentum 
cutoff imposed on the lattice. They were able to simulate 
a BEC collision with 150,000 bosons, using the positive- 
P representation for times long enough to obtain useful 
results, and thereby had an exact result to compare with 
the truncated Wigner method. The latter method pro- 
duced a "false halo" of particles in momentum space, de- 
pletion leading to unphysical negative densities beyond 
the halo, and accumulation of particles at low momenta 
- all in disagreement with the exact positive-P results. 
The Wigner method requires that initially empty modes 
of the system be represented by nonzero distributions, as 
if one half of a virtual particle occupied each mode. Evi- 
dently the truncated Wigner method treats these virtual 
particles as if they were real, in that a scattering event 
involving them can produce real populations of product 
modes. 

This is an ultraviolet divergence problem in that it be- 
comes worse as the momentum cutoff is increased. To 
obtain the most physically relevant results from a simu- 
lation, one must extrapolate to the continuum limit. It is 
in this limit, as the momentum cutoff approaches infinity, 
that the truncation errors are divergent. Clearly a full 
Wigner-Moyal treatment without truncation would not 
have these errors, but such a full theory with third-order 
derivatives also involves negative probabilities, which 
have no stochastic equivalent. 

We mention the projection method used with the trun- 
cated Wigner approach [l|| which amounts to an- 
other way to implement a cutoff, but does not solve this 
ultraviolet divergence problem. 

We mention here the projection methods as other tech- 
niques (not exact) for dealing with this problem . 

This discussion of problems with the truncated Wigner 
method is given as motivation for a Hybrid treatment. 
In future applications to multimode systems, we will in- 



vestigate whether the Hybrid method avoids these prob- 
lems. The large number of initial vacuum modes in a 
collision, for example, would be treated in the Hybrid 
method with the positive-P representation as phase-space 
variables set identically to zero. However this requires a 
detailed future investigation. The problem is absent in 
the pure positive-P method, although at long times very 
large sampling errors are found instead [9j. 



IV. THE POSITIVE-P METHOD 

The positive-P method involves an extension of the 
Glauber-Sudarshan P representation [13, HH from a sin- 
gle phase space to a doubled phase space, the same proce- 
dure that gives the doubled Wigner representation from 
the single. The defining equation (for a single-mode prob- 
lem) gives a representation of the density matrix in terms 
of a c-number function, P, and nondiagonal projection 
operators, Ap, both defined on a doubled phase space: 



d 2 a 



d 2 a + P(a,a + )A P (a,a + ), 



with 



A P (a,a+) 



|a)(a 



+*l 



(4.1) 



(4.2) 



Here |a) indicates a coherent state: a normalised eigen- 
state of the annihilation operator a. The effect of left- 
and right- multiplication of the density matrix by a and 



on P(a,a + ) can be deduced from Eqs. (|4. 1|4.2|) . The 



at 

proof involves an integration by parts in which bound- 
ary terms are assumed to vanish. The realm of validity 
of this assumption and the resulting effects on stochas- 
tic simulations are discussed at length by Gilchrist et al 
0. When the boundary terms vanish, the operator cor- 
respondences are: 



dp <-> aP(a, of 1 



pa 



a' p 



8a+ 



P(a,a A 



pa' 



a + P(a, a + ). 



(4.3) 



(4.4) 



(4.5) 



(4.6) 



As we noted in Table 1, all quartic Hamiltonian prob- 
lems, in the positive-P representation, give a true Fokker- 
Planck equation, with at most drift and diffusion terms: 

BP 1 

— = -d^A^a, a+)P) + -dMD^ia, a + )P), (4.7) 
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with 



and summation over fj,, v implied. Thus, for the many- 
boson Hamiltonian with two-body s— wave scattering 
terms [LT\ . no truncation of the positive-P equations is 
needed. 

The positive-P method solves two problems that oc- 
cur with single phase-space representations. First, if the 
distribution function is not guaranteed to remain real 
and non-negative, we cannot map the dynamics onto a 
stochastic simulation using standard methods. To deal 
with this, we use the feature of the representation that 
an infinity of different functions, P(a, a + ), may represent 
the same density matrix. We may choose, for the initial 
condition, a particular function 

P + (a,a+) = -L e ^l«-« + *l 2 (i(a + a +*)|p|i(a + a +*)), 

(4.9) 

that satisfies Eq . (|4. 1[) and is everywhere non-negative, 
as required. Alternatively, an initial pure coherent state, 
with p = I7X7I, can have a delta function representation, 
also positive: 

P(a,a + ) = 8 2 {a- 1 )8 2 {a + -7*). (4.10) 

The stochastic representation of this initial condition is 
simply 

a = 7, a + = 7*. (4-11) 

The second problem to deal with is that the diffu- 
sion matrix may not be positive semidefinite when writ- 
ten in the basis of real (x) and imaginary (y) parts 
(a x , a+, a y , a+). However, there is another symmetry in 
the positive-P representation, arising from the analyt- 
icity in a and a + of the nondiagonal projection opera- 
tor l|4.2p . that lets us make replacements to the real and 
imaginary parts of the derivative operators l|4.8p in (|4.3t - 
14.6ft . in just such a way that the resulting Fokker-Planck 
equation, relative to the component basis, has a positive 
semidefinite diffusion matrix With a positive initial 
condition and a true, positive semidefinite Fokker-Planck 
equation, the distribution is guaranteed to remain pos- 
itive. The standard method of mapping to a stochastic 
simulation also requires a diffusion matrix that is posi- 
tive semidefinite (all of its eigenvalues are non-negative) , 
so that stochastic equations are immediately derivable. 

The final step of mapping to stochastic differential 
equations involves first finding an A^-noise factorization 
of the diffusion matrix of the form 

JV 

D» v = ^B» n B vn . (4.12) 

This introduces another gauge degree of freedom that we 
will exploit later. Different choices of the factor matrix, 



B, that satisfy (|4. 12|) may provide stochastic simulations 
with widely different sampling error characteristics. 

The result of the adjustment of the diffusion matrix 
and this choice of the factor matrix is the set of Ito 
stochastic differential equations 

N 

da" = A"dt + B" n dw n , (4.13) 

n=l 

where the dw n are A" real Weiner increments [H] satisfying 
the stochastic average 

((dw n {t)dw m (t))) = S nm dt. (4.14) 

These SDEs, with appropriate initial conditions (equa- 
tions l|4.1ip for coherent states), are used to evolve a 
large ensemble of trajectories. The positive-P represen- 
tation is normally ordered, meaning that the most easily 
calculated quantum mechanical expectation values are of 
normally ordered operators. The formula for estimat- 
ing a normally ordered quantum mechanical expectation 
value as a stochastic average is: 

(tf m a n ) = ((a +m a n )). (4.15) 

We see from equation l|4.1ip that a coherent state can 
be represented initially with no noise in the positive-P 
representation. In this paper we will not embark on a 
detailed comparison of sampling error in the truncated 
Wigner, positive-P and Hybrid methods. However, we 
will take note of the number of trajectories needed, in 
each method, for an ensemble average to converge to a 
satisfactory result. All of our simulations were performed 
using xmds [12], and we used the built-in sampling error 
estimates of that program to judge convergence. 

The next part of our construction of the Hybrid repre- 
sentation involves writing the nondiagonal projectors for 
the positive-P representation in normally-ordered Gaus- 
sian form. The result of manipulating equation (|4.2p is 

A P (a, a+) =: g-C^ -«+)(«-«) . _ ( 4 16 ) 

V. PROBLEMS WITH THE POSITIVE-P 
METHOD 

A particular choice of the factor matrix, B, gives a 
set of stochastic differential equations (|4.13p that governs 
the evolution of the ensemble of trajectories. Unless this 
evolution is constrained in some way, trajectories may 
wander far from each other in phase space. Then the av- 
eraging over trajectories to estimate an expectation value 
may involve additions of many different, extremely large, 
numbers. No computer can calculate such an average 
without incurring a very large roundoff error. 

The result is the dramatic rise in sampling error that 
has been seen in some positive-P simulations. The 
growth in width of the distribution of trajectories often 
occurs over a short time scale, so that the sampling error 
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suddenly rises by many orders of magnitude, with the 
resulting growth of numerical errors. The simulation is 
of no value beyond this critical time. 

The problem can be caused by drift terms or noise 
terms, or a combination of both. A single-mode example 
to illustrate these problems is the anharmonic oscillator, 
with Hamiltonian 

1/ = wa'a + %a^a'flfl (5-1) 
and positive-P Ito stochastic differential equations: 

da = —i{u> + 2\Q. + a)adt + \J — 2ixdw (5-2) 

da + = + 2xa + a)a + dt + v / 2ix~dw + . (5.3) 

If we ignore the noise terms and choose a + a as a real 
number initially, the trajectory will be a circle in the 
complex a plane. But if a + a includes an imaginary part 
(from noise or from an initial condition other than the 
coherent state condition (|4.1ip ). either a or a + will spiral 
towards infinity, while the other spirals in towards the 
origin. The noise terms of the SDEs contribute to the 
problem, since they generally move a + a away from real 
values, thereby inducing the spiraling. 

Note that the single Wigner representation does not 
suffer from this problem because the real term \a\ 2 will 
always appear in place of a + a in the SDEs. 

Sampling error growth can be reduced or postponed 
by using our freedom to choose different factor matrices 
that give the same diffusion matrix, and by modifying 
the drift equations. Such methods are called stochas- 
tic gauge techniques [H, 0, HH]- However, while these 
are useful in single-mode examples, they are somewhat 
complicated when generalized to multi-mode cases. Also, 
we are interested in extending the time available for use- 
ful, error-free simulations to even longer time-scales than 
these methods can provide. 

A special case of problems with drift trajectories is 
when a trajectory is capable of reaching infinity in a finite 
time Q. This typically results in power-law tails in the 
distribution function, which violates the assumption that 
partial integration can be carried out. The simulation as 
it stands is then invalid beyond the singularity time. This 
problem can be dealt with using drift gauges [23]. Our 
examples will not fit into this category. 

We simulated the anharmonic oscillator in the positive- 
P representation to illustrate the sampling error problem. 
Figure 1 shows the X quadrature (X = -k(a + a')), with 
the choices ui — (for simplicity) and initial average num- 
ber N = 1. As we will do for all of our simulations, we 
plot results against a scaled time parameter, in this case 
Xt. (This is dimensionless when using H = 1). 

Deuar and Drummond have investigated various 
factors that affect the time for sampling error to be- 
come unmanageable in a multimode positive-P simu- 
lation. They have found that coarser spatial lattices, 




0.2 0.4 0.6 0.8 1 
Xt 

FIG. 1. X quadrature for single mode anharmonic oscillator 
vs xt'- positive-P method. Plotted are the ensemble average 
and the ensemble average ± sampling error estimate. The 
dashed line is the exact result. Parameters: w = 0, No = 1. 
Number of trajectories: 1,000. 

weaker interactions and lower particle densities all ex- 
tend the lifetime of the simulation. Of course the spatial 
lattice spacing can only be increased at the expense of 
systematic error, while the other two factors are fixed by 
the system being simulated. 

In the Hybrid scheme we will be using the positive-P 
representation only for the modes with lowest occupa- 
tions. Our test cases will investigate whether this delays 
the onset of large sampling error. 



VI. THE HYBRID METHOD 

The Hybrid method is designed to exploit a partic- 
ular feature of Bose-Einstein condensate systems: that 
a limited number of modes have very high occupation 
numbers. The method involves separating the physical 
system into modes that are, at least initially, highly oc- 
cupied (the condensed modes) and those that are lightly 
occupied (the output of an atom laser or the products of 
a BEC collision). Then we intend to use different rep- 
resentations to treat different modes, treating the highly 
occupied modes with a form of the Wigner representa- 
tion and the lightly occupied modes with the positive-P 
representation. 

Use of the Wigner representation for the highly occu- 
pied modes will in general simplify the structure of the 
resulting diffusion matrix. In a simple two-mode model 
discussed in Section VII, we will see that this allows us 
to delay the rapid growth of sampling error. And by 
not using the Wigner representation for the potentially 
very large number of lightly occupied modes, we intend 
to avoid the false halo problem. 
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Our first task is to show that we can consistently 
use two different representations on different modes, cor- 
rectly describing interactions that couple these different 
modes. For general interactions of this sort, diffusion 
terms involving the Wigner modes are inevitable. To be 
able to construct a positive semidefinite diffusion matrix 
in the general case, we will have to use a doubled phase 
space throughout. 

Now we can exploit the similarities in equations (|2,6p 
and (|4. 16|) to define a Hybrid representation with a par- 
ticularly simple notation. We suppose that a system has 
modes labelled m = 1, . . . , M . These modes are to be 
treated with the Wigner representation or the positive- 
P representation depending on whether a parameter r m 
takes the value: 

r m = 1 for positive-P, (6-1) 

r m = 2 for Wigner. (6-2) 

The nondiagonal projection operator is a direct product 
of terms for each mode, m: 

M 

k H {a;r) = JJ r m : e ~r m & -<*+)(&-*) . _ ^ 

m— 1 

Then the Hybrid representation of the density matrix 
becomes 

p = J d 4M dP H (d;r)A H (a;r), (6.4) 

where a — (ai, af , . . . , aw, a^) and f— (r\, . . . , ru )• 

Note that the use of the parameter r for these doubled 
phase space representations is very much like Glauber 
and CahilPs use of the parameter s to span antinor- 
mally ordered, symmetrically ordered and normally or- 
dered single phase space representations. The connection 
between the two schemes follows by taking 

s — 1 . 
r = (6.5) 

s 

and mapping doubled to single phase spaces. 

The applicability of the positive-P method depends on 
the two results that we mentioned in Section IV. First, 
for any initial density matrix, it is possible to choose 
a phase-space distribution that is everywhere real and 
non-negative (using equation (|4.9p ). Second, it is always 
possible to cast any diffusion matrix into a form that is 
equivalent with respect to physical predictions and that 
is positive semi-definite in the basis of real and imagi- 
nary parts of the phase-space variables. Corresponding 
results must hold for any Hybrid representation in order 
for that method to be usable. We were able to prove both 
assertions. 

First, we found an integral transform that takes a 
positive-P distribution to a doubled Wigner distribu- 
tion representing the same density matrix. We show the 



single-mode case: 

W(a,a + ) = ± J e-^-^'? P{\W+X),\W-X)), 

(6.6) 

with 

t/j = a + a + *, x — a — a^". (6-7) 

Note that P has four independent real parameters, but 
the integration is over only two degrees of freedom. The 
single-mode case is shown but the extension to the multi- 
mode case is straightforward. Since the kernel is positive, 
the transform can be used to take the initial positive dis- 
tribution l|4.9p to an everywhere positive doubled Wigner 
distribution. Extension to the case with many modes 
treated by different representations proves the first as- 
sertion. 

The proof of the second assertion is exactly like the 
textbook proof for the positive-P representation, since 
the derivative equivalences 



d 

da 


d 

<— > <— > 

da x 


d 
da y ' 


d 
da+ 


d 

<— > — — <— > 

dat 


d 
dat 



are the same as their positive-P counterparts. 

Use of the Hybrid method is simple for few- mode prob- 
lems. For the mapping of the evolution equation l|3.fp 
for p to a Fokker-Planck equation, we use either the 
Wigner (|2.7H2.10p or positive-P H4.3M.6p operator cor- 
respondences as appropriate for each mode. In gen- 
eral there will be terms with three derivative opera- 
tors for quartic Hamiltonians. (Terms with four deriva- 
tives always cancel.) For each application, we must de- 
cide whether truncation of these terms, to produce a 
drift /diffusion problem, is valid. Scaling arguments like 
those applied to the Wigner method l|3.3p can be used 
here. In problems involving both highly occupied modes 
and lightly occupied modes, there may occur problematic 
three-derivative terms from mutual interaction of those 
modes. 

A feature peculiar to the Hybrid method is that there 
will appear what we call interface noise: there will be 
diffusion terms that are proportional to the difference of 
r values for different modes, that would vanish if those 
modes were treated with the same representation. 

The mapping to stochastic differential equations uses 
the same rule as is used for the positive-P representation: 
if a generally complex matrix B provides a factorization 
D = BB T of the diffusion matrix, then the Ito stochastic 
differential equations can be chosen as 

N 

da^ = A^dt + ^2 B^ n dw n , (6.10) 
n=l 

where fi labels the components of the vector of phase 
space variables a = (a%, a^ , . . . , au-, &m) an d the dw n 
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are N real, independent Weiner increments. We note 
that the freedom of choice of a factor matrix, B, intro- 
duces a gauge degree of freedom that may allow us to 
reduce sampling error in simulations. 

The relation between physical expectation values and 
stochastic averages will take new forms in the Hybrid 
representation. Here an observable may be a product 
of factors to be treated with the symmetrically ordered 
Wigner representation and others to be treated with the 
normally ordered positive-P representation. So, for ex- 
ample, in Section IX we will need to calculate an expec- 
tation value as 

(N a Y b ) = (Iat*(5 - 6t)) = {{^ + a - \){P ~ /3 f ))), 

(6.11) 

where the a mode is treated with the Wigner represen- 
tation while the b mode is treated with the positive-P 
representation. 



VII. TEST CASE: COUPLED ANHARMONIC 
OSCILLATORS 

As a first test of the Hybrid method, we simulated the 
behavior of two coupled anharmonic oscillators, with a 
coupling that preserves the individual mode occupations. 
The Hamiltonian is 

H = u a a' a+Xa& & a&+&bb b+Xbb b 'bb+ga 'a& (7.1) 

We used an initial coherent state for the a mode (with 
high mean occupation iV a o = 100) and for the b mode 
(low mean occupation Nbo = 0.01). We set uj a = ^>b = 
for convenience and used Xa = Xb = 9 = lj which sets 
the scale for the time variable. 

This model is meant to resemble just a few terms of 
the much larger multimode Hamiltonian for a Bose gas 
with s-wave scattering terms. 

Note that we have chosen a model system in which the 
a occupation remains constantly large, while the b occu- 
pation stays small. The Hybrid method can be used in 
cases where these numbers are not conserved, and gives 
good results when the occupations of the modes remain 
high and low over the interaction time, respectively. Re- 
sults from this category will be presented in a later work. 

We simulated this system using the Hybrid method 
and, for comparison, the truncated Wigner method and 
the positive-P method. We were also able to obtain an 
exact solution for coherent state initial conditions, as did 
Chaturvedi and Srinivasan (26| . 

We insert the Hybrid representation (|6.4p of the den- 
sity matrix into the evolution equation (|3.1|l . An inte- 
gration by parts, justified in this case, amounts to using 
the Hybrid operator correspondences (|2.7H2.10| 14.314.6(1 . 
This gives an equation of the form 



MP, 



d a a^^-K H {a,r) = I <PaC(a, r)P H (a, r)A H (a, r), 



Ot 



(7.2) 



where £ is a linear, differential operator that acts on P#. 

We extract a Fokker-Planck equation for Pjj from 
(|7.2p . keeping all terms, including third-order derivative 
terms. (We note that for doubled phase-space represen- 
tations this choice is not unique.). We find 



.dP H 



d_ 

da 
d 



{2 X a(a + a-l)+gl3 + l3}aP H 



da 
d_ 

dp 




-{2 Xa {a + a- 1) 



gp+P}a+P H 
1. 



{2 X bP + P + g(a + a~-)}pP H 



1. 



u i+ {2 X bP + P + g{a + a~-)}p + P H 



d 2 32 d 2 

+XbgpP PH-Xbg^ 



P +2 P H 



9 d d g d d + 

'2d^dj3 aPPH+ 2d^W aP PH 



g d d + g d 

2 da+ dp 



d 



9_d ^_(^_ g 

Ada da+ l dp P ' 

Xa r d 2 d 



dp 

d 2 



2 da+ dpi 

\p + }P H 



■a+P+P H 
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2 da 2 da + da +2 da 



a+}P H - (7.3) 



We use the convention that the derivative operators act 
on all factors to the right. 

To apply the scaling argument discussed in Section III, 
we write the third-order derivative terms above (which we 
call T 3 ) in terms of the scaled phase-space variables 



u = a/ 




^a+/ 




(7.4) 
(7.5) 



Then T3 becomes 

1 g d 



T, = 



8 



N aa 4 du du+ 

1 Ya r d 2 



{- 



I 9 

dv 

d 



-v+}P H 



N a0 2 l du 2 du^ 



dv 

d 2 d 

— 1 

du+ 2 du 



■}Ph.(7.6) 



We expect, in the stochastic simulation of this problem, 
that there will be a finite time scale over which a and a + 
will remain distributed close to order y/N a o in magnitude, 
while P and P + remain near \fNm- Our first simulation 
will stay within this time region. Over that time scale, 
the third-order derivative terms will make a negligible 
contribution to Pr compared to the drift terms (first- 
order terms which scale like A^o) and the diffusion terms 
(second-order terms which scale like 1). 

After we truncate these terms, the Fokker-Planck equa- 
tion has drift vector (in the basis (a, a + ,P,P + j) 



A 



( -i{2xa(a + a-l)+gP+P}a \ 
+i{2xa(a + a-l)+gP+P}a+ 
-i{2 X bP + P + g{a+a-\)}P 
\+i{2 Xb p + P + g{a+a-\)}p+ 



(7.7) 
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and diffusion matrix: 



D = 






-fa/3 
-fa/3+ 






-f £*+/?+ 



-fa/3 




-f a+3+ 



(7., 



Because of the use of the Wigner representation for 
the a mode, this diffusion matrix differs from the one 
resulting from a pure positive-P treatment in the absence 
of terms — 2ix a a 2 and +2ix a a +2 in the first two diagonal 
spaces, respectively. 

We were able to construct a factorization of the dif- 
fusion matrix (|7.8p by first treating the diagonal terms 
and then recognizing a simple structure in the remaining 
matrix. The following factor matrix requires only four 
real noises in the SDEs: 



D 





(o 




























b 


i/3 






















v) 






(0 





a 


ia \ 




-ig 








—a 


+ —ia + 










P 














(7.9) 



The resulting SDEs produced the results shown in Fig- 
ures 2 and 3. We calculated the expectation values of the 
quadrature operators X a — ^{a + a)) and Xb = ^(b+w). 
The simulation was clearly stable over the time scale 
shown and gave results in excellent agreement with the 
exact solution. We will refer to the method used here as a 
gauge Hybrid method, since it relies on being able to find 
a diffusion gauge (a useful factorization of the diffusion 
matrix) . 



x" o 




PIG. 2. X quadrature for mode a vs x a t: coupled 
anharmonic oscillators treated with the gauge Hybrid 
method. Plotted are the ensemble average, the ensemble 
average ± sampling error estimate and the exact solution. 
Parameters: u a = cut = 0, \a = Xb = g, N a o = 100, 
iYi,o=0.01. Number of trajectories: 10,000. 



0.05 








-0.05 



FIG. 3. X quadrature for mode 6 vs \ a t: coupled 
anharmonic oscillators treated with the gauge Hybrid 
method. Plotted are the ensemble average, the ensemble 
average ± sampling error estimate and the exact solution. 
Parameters: u a = oJb = 0, % a = x& = 9 = 1, N a0 = 100, 
N b0 = 0.01. Number of trajectories: 10,000. 



When we simulated this same problem using the trun- 
cated Wigner method, the results were nearly indistin- 
guishable from Figures 2 and 3 (using 150,000 trajecto- 
ries), so we do not display them here. With regard to 
this first test, we have not yet established superiority of 
the Hybrid method over the truncated Wigner, except 
to note that the Hybrid method requires far fewer tra- 
jectories to attain a given accuracy. In Section VIII we 
will explore a different region of parameter space and in 
Section X we will calculate a higher-order moment in the 
same system. In both cases, we will see results that show 
a clear distinction between the methods. 

We also simulated this problem with the positive-P 
method. Sampling error rose to very large values at about 
t = 0.04, after just one oscillation of the quadratures. 



Analysis of the third-order derivative terms from the 
pure Wigner calculation, similar to the above analysis for 
the Hybrid method, shows terms that scale as l/iV^o and 
so cannot be justifiably neglected. 

The mechanism at work in stabilizing the Hybrid sim- 
ulation over limited times is as follows. With this choice 
of gauge, the stochastic differential equations keep the 
quantity a + a fixed, for each trajectory, at its initial 
value. These values, selected by the stochastic Wigner 
initial condition of the form l|2.11l I2.12[) , will always be 
real and close to iV a o. The quantity /3 + (3 starts at 7V{,o 
then acquires an imaginary part, but its magnitude is 
kept of order Nbo over the simulation time. 

Further inspection shows that the magnitudes of a and 
a + will remain near \/N a o while those of (3 and /3 + remain 
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of the order of V^bo over the simulation time. (These 
estimates were used to justify our neglect of the third- 
order derivative terms in equations (|7.4l I7.5( 17. 6§ .) So 
the drift terms are dominated by the factors of a + a and 
spiraling is negligible. 

Over a short time At, the relative sizes of the drift 
and diffusion increments, for z one of the phase space 
variables, are given by 

Drift: Az~N a zAt 



the fact that the quadratures are strongly damped be- 
fore the neglect of terms (for both methods) and sam- 
pling error growth (for the Hybrid method) can become 
important. We lowered the mutual interaction strength 
between the two modes, relative to \a — Xb> by set- 
ting g/xa = 0.0001. This greatly extended the damping 
time for the b mode, allowing us to see differences in the 
predictions of the gauge Hybrid and truncated Wigner 
methods. 

The results are shown in Figure 5. 



Diffusion: Az ~ z \/~Ab 

(with Xa = Xb — 9 — 1)- So diffusion is, in this exam- 
ple, negligible compared to drift over the time scale of 
interest. 

We calculated the quadrature X a in our model to 
longer times, with results shown in Figure 4. (To obtain 
the qualitative features rapidly, we used, in each 
lower number of trajectories than we used in our previous 
simulations.) The exact result showed a recurrence cen- 
tered on t — it. The gauge Hybrid method showed large 
sampling error before that time, starting at about t = 2.5. 
The truncated Wigner method was also unable to predict 
this recurrence, showing instead a quadrature remaining 
close to zero. Recalling that the pure positive-P treat- 
ment suffered large sampling error after about t = 0.04, 
we see that the gauge Hybrid method extended the useful 
simulation time by a factor of 60. 




FIG. 4. X quadrature for mode a calculated to longer 
(dimensionless) times \ a t. (a) Exact result, (b) Gauge 
Hybrid result: 100 trajectories, (c) Truncated Wigner result: 
10,000 trajectories. Parameters: oj a = <^b = 0, \a — Xb — 9, 
N a0 = 100, Nw, = 0.01. 



VIII. WEAK COUPLING 

In the previous example, both the gauge Hybrid 
method and the truncated Wigner method are aided by 
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FIG. 5. Comparison of the gauge Hybrid and truncated 
Wigner methods in predicting the X quadrature for the b 
mode of the coupled anharmonic oscillators, with weak 
coupling. Quadratures plotted against XoX- Gauge Hybrid: 
10,000 trajectories. Truncated Wigner: 15,000 trajectories. 
Parameters: u a = oj b = 0, \a = Xb, g/Xa = 0.0001, 
N a0 = 100, Nw = 0.01. 



We see that the truncated Wigner method fails from 
Xat = 0, consistent with our expectations for a system 
with a very lightly occupied mode. The gauge Hybrid 
method performs well until about Xat — 1.5, when it is 
overwhelmed by sampling error. 

IX. FURTHER TRUNCATION 

In the examples we have seen so far, truncation of 
terms in the Hybrid method has not prevented it from 
attaining excellent agreement with the exact solutions at 
early times, even when dealing with a very lightly oc- 
cupied mode. The method is, however, clearly limited 
by the growth of sampling error. In this section, we try 
a simple adjustment to the equations to try to extend 
useful results to longer times. 

When phase space distributions grow wide in uncon- 
strained directions, the trajectories sampling those distri- 
butions are widely spread and the calculation of expec- 
tation values becomes a great numerical difficulty. To 
understand the meaning of the widths in "unconstrained 
directions," we note that we could estimate the spread of 
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our Hybrid distribution by calculating all the stochastic 
averages ((aW)), where a 1 is a real or imaginary part 
of a (defined after equation (|6.10p ). Some of the linear 
combinations of these averages, such as ((a + a>)), are con- 
strained to approach physical predictions as the number 
of trajectories grows large. Widening of the distribu- 
tion in the other directions will increase sampling error, 
but may be reduced using the gauge freedoms of doubled 
phase-space representations, or other methods. 

Spiraling of the drift trajectories is one source of 
spreading that we have identified, and that we have al- 
ready partially controlled using our choice of gauge. For 
our gauge SDEs, the quantity a + a remains completely 
real for all times, and thus does not cause spiraling in the 
drift equations. Not so the quantity /3 + /3, which starts 
with a purely real value but can immediately develop an 
imaginary part from the influence of the noise terms. 

We tried a further truncation of our gauge Hybrid 
equations, making the replacement 



(9.1) 



In future applications, if a + a is not constrained, we pro- 
pose to also try the truncation 



a + a —>■ Ke(a + a). 



(9.2) 



We saw good short-time behavior from this truncated 
Hybrid method, equaling that of all the other methods. 
At longer times the method was unable to predict the 
recurrence, showing quadratures staying close to zero. 
But the sampling error remained at a manageable level to 
t = 5.0. In future work, we will investigated whether this 
somewhat ad hoc truncation can be used as a simple way 
to extend simulations to longer times without incurring 
excessive systematic error. 



X. TEST CASE: QUANTUM NONDEMOLITION 
MEASUREMENT 



Our first test case showed the Hybrid method — with 
a diffusion gauge choice and with a further truncation — 
able to successfully simulate an interacting system be- 
yond the time at which the positive-P method became 
unusable. But the Wigner method was able to give 
equally good results on the same system. (A distinction 
was found in the weak coupling case.) Here we investi- 
gate a different observable — a higher order moment — in 
the same system, and find the results more sensitive to 
the choice of method. 

The concept of quantum nondemolition measurements 
[13, arose from the need for a way to measure the very 
small displacements of a gravitational wave detector that 
are expected to occur from the passage of a gravitational 
wave. Repeated measurements of position, to high ac- 
curacy, would be required to distinguish the signal from 



other effects. Quantum mechanics sets limits on schemes 
to measure those small displacements. Measurement of a 
position observable with a finite uncertainty may produce 
a state in which the uncertainty in position grows after 
the measurement. At later times, when another measure- 
ment of position is performed, the uncertainty would be 
larger that the desired maximum. 

Instead, measurement of a conserved observable, such 
as the momentum of a free particle, can be repeated an 
arbitrary number of times without causing the uncer- 
tainty to increase. The quantum nondemolition (QND) 
measurement scheme involves choosing an appropriate 
conserved observable (in a probe beam) that can give in- 
formation about the signal of interest after the signal and 
probe interact. 

A QND scheme can be constructed from our model 
of interacting anharmonic oscillators f29j. We suppose 
that the bosons in question are now photons, and that 
they can interact with each other in a suitable nonlin- 
ear medium, such that our number-conserving interac- 
tion Hamiltonian gives a toy model of the dynamics. Of 
course a fuller description of the dynamics would involve 
propagation in space, dispersion and other factors (3(J. A 
lightly occupied signal beam and a highly occupied probe 
beam interact in the medium. Phase information will be 
exchanged between them, while their individual number 
distributions are conserved. 

In one QND scheme, the conserved QND observable is 
taken as the photon number, N a , in the highly occupied 
probe beam. The signal is the phase quadrature of the 
lightly occupied beam, Yj, = — ^(w — b). 

We suppose that the interaction between signal and 
probe lasts only for a short time, as would be the case for 
two short pulses interacting in an optical fiber. We make 
the interaction cease when the magnitude of the corre- 
lation function reaches its first maximum. This means 
that we use the previous Hamiltonian of equation (|7.1|1 , 
except with g = 1 for t < to, and g = for t > t , where 
to = 0.1 in our example. 

We calculate the correlation function between probe 
and signal, a measure of the potential success of the mea- 
surement scheme: 



C(N a ,Y b ) 



(N a Y b )-(N a )(Y b ) 
vHN a )Vi{Y h ) 



where 



(10.1) 



(10.2) 



is the variance of an operator ft. We set u> a = for conve- 
nience, since this will remove a high frequency variation 
from our expectation value. Likewise, we set uj b = —N a ag 
to obtain a slowly varying expectation value. This latter 
choice is equivalent to a particular choice of local oscilla- 
tor frequency in the homodyne detection of Y b . 

Figure 6 shows the correlation function calculated with 
two different phase-space methods and compared to the 
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exact result. The gauge Hybrid method shows excellent 
agreement with the exact result. In contrast, the trun- 
cation of the Wigner method evidently removes terms 
that are needed to correctly predict the correlation func- 
tion at times after the interaction ceases. The tendency 
of the truncated Wigner method to give worse results 
when predicting higher-order moments was investigated 
by Drummond et al 
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FIG. 6. Comparison of methods for determining the 
correlation between N a and Yb for a QND scheme, vs \ a t. 
Results are shown for the Hybrid method with a diffusion 
gauge (50,000 trajectories) and the truncated Wigner 
method (50,000 trajectories), compared to the exact result. 
Parameters: N a0 = 100, Nto = 0.01, uj a = 0, w;, = —N a og, 
Xa = Xb, g/Xa = 1 for Xat < 0.1, g = for X at > 0.1. 



of rapid sampling error growth by a factor of 60 com- 
pared to the positive-P method. Further investigations 
will focus on many-mode systems to see whether these 
advantages over the earlier methods can be maintained. 

It is interesting to note here that our results show that 
a very natural application of the Hybrid method is to 
systems of two different types of particle with interactions 
that conserve individual species numbers. This presents 
a natural framework to investigate quantum Brownian 
motion, which will be treated in subsequent work. 
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APPENDIX: EXACT SOLUTIONS 



The Hamiltonian 

H = uj a a)a + XaO^a^aa + uj b wb + Xb&b'bb 



gci'aPb, 



(A-l) 

describing two coupled anharmonic oscillators, can be 
written just in terms of the number operators, N a = a) a 
and N b = b'b, as 

H = u a N a +u b N b + Xa {N 2 a ~N a ) (A-2) 
+Xb{N h 2 - N b ) + gN a N b . 

So the number states 



XI. CONCLUSIONS 

We have shown that, for a stochastic phase-space treat- 
ment of a multimode system, it is possible to use the 
doubled Wigner representation for some modes and the 
positive-P representation for the remainder. We tested 
our method on a system of two coupled anharmonic oscil- 
lators, one with a mean occupation that remained at 100, 
the other with a mean occupation of 0.01. The method 
was able to simulate the evolution of quadrature expec- 
tation values for times far beyond where the positive-P 
method suffers a rapid growth of sampling error. Results 
were in excellent agreement with the exact solution. 

While the truncated Wigner method performed as well 
as the Hybrid method when calculating these quadra- 
ture observables (over a finite time), for the calculation 
of a higher order moment corresponding to a QND ex- 
periment there was a very clear advantage of the Hybrid 
over the truncated Wigner. The latter results contained 
a large systematic error, while the Hybrid result was in 
excellent agreement with the exact result. 

At least as applied to this system with a small number 
of modes, the Hybrid method was able to delay the onset 



\n a n b ) 



rIO) 



(A-3) 



V" Q ! V n b 

are eigenvectors of the Hamiltonian with eigenvalues 

E(n a ,n b ) = u a n a + u b n b + Xa(nl - n a ) 

+Xb{nl - n b ) + gn a n b . (A-4) 

We consider an initial state that is a coherent super- 
position of the number states (|A-3[) of the form 



I7a7fc; 



n 00 Tt-b 

n a =l VW„^ V^fc! 



E 



(A-5) 

where 7„ = y/N a o, jb = y/N b o and iV a0 and N b0 are the 
average occupations of the modes. 

We are interested in observables, fl, that are simple 
combinations of a small number of creation and/or anni- 
hilation operators. These have simple matrix elements, 
(n' a n b \Cl\n a rib) , between the number eigenvectors. All 
terms will be proportional to Kronecker deltas of the form 
S n ' a ,n a +rn a S n ' b ,n b +m b , for to q and rn fc integers. Then the 
expectation value of such an operator in the state vector 
produced by time evolution of (|A-5|1 will always reduce 



13 



to a double sum (over n a and rib), with the unitary time 
evolution producing a known phase factor inside the sum. 

The time-dependent expectation values of the quadra- 
ture operators X a = |(a + aJ), Y a — ^(a — a 1 ), X b = 
\{b + tf) and F b = ^(b - for the initial state 1|A-5|1 . 
can then be evaluated as sums over n a and n b . We find 
the results 

(X a (t)) = y^ e -{^o(l-co S 2 Xat )+Ar bo (l-cosc,t)} ( A _g) 

x cos{uj a t + N a0 sin 2\J + N b0 sin gt}, 

(Y a (t)) = — ^/ N a Qe^^ Na0< - 1 ^ cos2XatS>+Nbol - 1 ^ cos 9t ^ 

x sin{cj a t + N a0 sin2x a i + N b0 sin gt}. (A-7) 

For (Xb(t)) and (Yb(t)), we make the replacement a «-> 6 
in expressions (| A-6(l and (j A-T[) . respectively. 

Our model of a QND measurement has the feature 
that the coupling strength, g, is constant up to a time, 
r, and vanishes after that. This is meant to model 
two light pulses that cease to interact after they no 
longer overlap within an optical fiber. The evolution op- 
erator for the resulting time-dependent Hamiltonian is 
U (t) = exp — J Q H(t')dt' . For times up to r, we evaluate 
the following: 



V(Y b ) = (f b 2 ) - (Y b f 



(N a Y b (t)) = -N a y/N~ b e- NMl - COS9t) 



xe 



-N b0 (l-cos2 X bt) 



(A-8) 



x sin{(w b + g)t + N a0 sin gt + N b0 sin2x 6 i}, 



_ 1 Nboe ~{Na(i-cos2gt)+N b (i-cosi Xb t)} x 

x cos{2(w fc + Xb)t + N a sin 2gt + N b sin 4x&i} 

_ Nbe -2{N a {l~cosgt)+N b (l-cos2 X bt)} x ^_g% 

x sin 2 {ui b t + N a sin gt + N b sin 2xbt} 
+ 7 + l N bo- 



For times beyond r, we make the replacement gt —* gr in 
(|A-8j) and (|A-9[) . This rule applies also to the expectation 
values 1|A-6|) and (|A-7|) for t > t in this QND scheme. 

Finally we note that the individual particle numbers 
are conserved under this Hamiltonian, and the coherent 
state initial conditions (| A-5[) give 



(N a ) = iV a0 , (N b ) = N t 



M) • 



(A-10) 



V(N a ) = N a0 , V(N b ) = N, 



bO- 



(A-ll) 
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